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ABSTRACT 

We calculate the radial diffusion coefficient for a passive contaminant in an accretion 
disc which is turbulent due to the action of the magnetorotational instability. Nu- 
merical MHD simulations are used to follow the evolution of a local patch of the disc 
using the shearing box formalism. A separate continuity equation for the mass fraction 
of contaminant is integrated along with the MHD sytem, and radial profiles of this 
fraction are obtained as a function of time. Solutions of a linear diffusion equation 
are fitted to the numerical measured profiles of the contaminant, treating the diffu- 
sion coefficient D as the fitting parameter. At early times, the value of D is found to 
vary, however once the contaminant is spread over scales comparable to the box size, 
it saturates at a steady value. The ratio of D to the transport coefficient of angular 
momentum due to shear stress is small. If D can be used as a proxy for the turbulent 
magnetic diffusivity, the effective magnetic Prandtl number Peff = i^/D (where ly is 
the coefficient of "effective viscosity" associated with shear stress) would be large. 

Key words: accretion: accretion discs - magnetohydrodynamics 



1 INTRODUCTION 

In order to understand the morphology and chemical com- 
position of protoplanetary accretion discs, the radial mixing 
of vapour and solid material in accretion flows needs to be 
characterised (Morfill 1983) . The distribution of dust across 
a turbulent disc, as a consequence of both advective and dif- 
fusive modes, will determine the mechanism which is most 
favourable for the assemblage of particles into planetesimals, 
meteoritic bodies and eventually planets. For example, there 
is evidence for radial upstream movement of nebular gas 
from the coexistence of calcium-aluminum inclusions and 
chondrules in meteorites (Boss 1998). Studies by Takeuchi 
& Lin (2002) and by Keller & Gail (2004) indicate that, for 
vertically isothermal laminar discs, there may exist an equa- 
torial outflow capable of transporting material from the hot, 
inner regions to the cold, outer regions, thereby affecting the 
composition of bodies located in the latter zone. The mate- 
rial is assumed to be well coupled to the gas, which is true 
for particles whose size is smaller than the mean free path 
of the gas molecules. However, this picture may be changed 
if angular momentum transport in the disc is mediated by 
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magnetohydrodynamical turbulence driven by the magneto- 
rotational instability (MRI) (Balbus & Hawley 1991). 

The transport of a passive scalar has been used as a 
criterion to study diffusive processes in turbulent flows (e.g. 
Brandenburg, Kapyla & Mohammed 2004 and references 
therein). Research on the behaviour of passive tracers in 
turbulence has shown that there is a strong variability of dis- 
sipation and mixing rates of both scalar and velocity fields, 
due to small scale intermittency. It has also been shown 
that the central assumption of Kolmogorov scaling, that at 
small scales the scalar and velocity fields are isotropic in the 
limit of infinite Reynolds and Peclet numbers, breaks down 
(Warhaft 2000; Sreenivasan 1991). Although a fascinating 
subject, in this work we will not address many important 
issues that belong to the vast field of passive scalars. 

Clarke & Pringle (1988) analysed the evolution of the 
concentration of a trace contaminant in a steady, Keple- 
rian disc. They point out a subtlety related to the identi- 
fication of a diffusion coefficient with a viscosity for tur- 
bulent accretion. Despite the fact that both are associ- 
ated with velocity fluctuations in the turbulence, their ef- 
fects are distinct: diffusion tends to produce a uniform mix- 
ing of material, whereas viscosity does not yield a uniform 
distribution of angular momentum. They found that the 
contaminant behaviour depends on the ratio of the diffu- 
sion coefficient to viscosity; only when this ratio is near 
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or greater than unity docs diffusion take place apprecia- 
bly. Heyvaerts ct al. (1996) have shown that in an accre- 
tion disc in which the effective viscosity is due to turbulence 
generated by the Balbus-Hawley instability, the associated 
magnetic Reynolds number is small, whereas the magnetic 
Prandtl number Pm = v/r) (where t] and u are the mag- 
netic resistivity and the effective viscosity, respectively) is 
not necessarily so. 

In this context, we wish to investigate the role that tur- 
bulent diffusion plays locally in the transport of a passive 
tracer, within a small portion of a weaJily magnetised ac- 
cretion disc in which the MRI operates. We use numerical 
MHD simulations which adopt the so-called "shearing box" 
model, discussed in detail by Hawley, Gammie & Balbus 
(1995, hereafter HGB), a useful scheme for the study of lo- 
cal processes in accretion flows. We extend the numerical 
scheme to explicitely follow the time evolution of a passive 
contaminant in the flow to determine the rate of turbulent 
mixing. We use the ideal MHD approximation to simplify 
the problem; this requires the magnetic field be well coupled 
to the gas. Although the structure of protoplanetary discs is 
not well known, it is likely that far from the central star, the 
central regions of the disc arc dense and cool enough that 
non-ideal MHD effects such as Ohmic dissipation, ambipolar 
diffusion, and the Hall effect will significantly alter the MHD 
(Blaes & Balbus 1994; Wardle 1999; Sano et al. 2000; Stone 
et al. 2000; Balbus & Terquem 2001; Salmeron & Wardle 
2003). Recent studies of the MRI in non- ideal MHD indi- 
cate turbulence and transport my be strongly suppressed if 
the magnetic Reynolds number is close to one (Hawley & 
Stone 1998; Sano & Miyama 1999; Fleming, Stone, & Haw- 
ley 2000; Sano & Stone 2002). It is possible that even ff 
the MRI is suppressed in the central regions of the disc, the 
surface layers may be ionized by non-thermal effects and re- 
main turbulent (Gammie 1996; Glassgold et al. 1997; Flem- 
ing & Stone 2003). Of course, by adopting the ideal MHD 
approximation, this study cannot be applied to regions of 
the disc where the ionisation fraction is so low that it is sta- 
ble to the MRI: the turbulent transport of contaminants in 
these regions will be essentially zero. In regions where the 
MRI is suppressed but not eliminated, the saturation am- 
plitude of the turbulence will be reduced, thereby reducing 
both the coefficient of effective viscosity v and the turbulent 
coefficient D. Thus, the ratio of the two (the effective mag- 
netic Prandtl number Pcb) may be independent of the level 
of saturation of the instability. Studying diffusion in realis- 
tic, global models of protoplanetery discs including non-ideal 
MHD and the vertical and radial dependence of the ionisa- 
tion fraction is a challenging problem and beyond the scope 
of this study. 

Previous work has investigated diffusive processes in dif- 
ferent magnetohydrodynamic contexts. Yousef et al. (2003) 
performed forced turbulence simulations in which they ob- 
tain a value of the turbulent magnetic Prandtl number of 
about unity, regardless of the value of the microscopic mag- 
netic Prandtl rmmbcr. Brandenburg et al. (2004) study non- 
local aspects of turbulent transport by adding a second- 
order time derivative to the equation for the concentration 
of a passive contaminant, namely Pick's law, which states 
the proportionality of the flux of passive scalar to the neg- 
ative concentration gradient. This extra term results in a 
damped wave equation for the concentration. Our investi- 



gation is similar, except turbulence in our models is driven 
by a linear instability in the background flow, rather than 
by an assumed forcing. Moreover, the orbital dynamics in an 
accretion flow strongly affect the velocity and magnetic field 
correlations in the turbulence, which may directly affect the 
turbulent diffusion. 

This paper is organised as follows. In Section 2 the nu- 
merical method employed is described, including the method 
by which the turbulent diffusion coefficient is calculated. Re- 
sults of a variety of simulations arc described in Section 3, 
while a discussion and conclusions are presented in Section 
4. 



2 METHOD 

2.1 Numerical Method 

We use a three-dimensional version of the ZEUS code (Stone 
& Norman 1992a;b) in which the shearing box model of 
HGB is implemented. Briefly, a coordinate system is cen- 
tered at a fiducial radius ro in a small patch of a disc, which 
has linear dimensions much less than ro. This new system 
corrotates at an angular velocity fio = f2(r'o). The x,y,z 
coordinates correspond to cilindrical r, ip, z coordinates, re- 
spectively. The resulting MHD equations then incorporate 
the effects of Coriolis forces and tidal stresses (see below). 
The numerical implementation of the shearing box involves 
boundary conditions such that the computational domain is 
strictly periodic in the y, z directions at all times, whereas 
in the x (radial) direction it is quasiperiodic: when a fluid 
element exits at one radial boundary, it reappears at the op- 
posite radial boundary with a position and a velocity given 
by the shearing motion. 

The passive contaminant is introduced in the code as a 
fraction of the gas density. In addition to the MHD equations 
for the shearing box. 



dp 
dt 
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where the symbols have their usual meanings, the code also 
solves a continuity equation for the contaminant fraction: 

dfp 



+ V • (/pv) = 



(6) 



at 

where / is a number between and 1. 

As in HGB, we use the standard values P = Pq — 10~® 
(initial pressure), H = 1 (box height), and Q, — 10^"^. Our 
runs adopt the same parameter values as those of model Z4 
of HGB (in particular, we use a uniform strength vertical 
field with plasma parameter /3 = 400), except that we use 
a resolution of 88 x 184 x 88 grid zones, almost three times 
greater in each dimension than that used in Z4. The box size 
is 1 X 27r x 1. 



Diffusion coefficient of a passive contaminant in a local MHD model of a turbulent accretion disc 3 



2.2 Simulations 

HGB showed that turbulence develops after a few orbits 
in the evolution of the shearing box (one orbit=27r/f2). We 
first run the code for 16 orbits, to ensure that turbulence 
has set in and is able to stir the contaminant efficiently. 
At t = 16 orbits the code is restarted and the contaminant 
is introduced; we set / = 1 for —0.15 < x < 0.15 and 
for all y and z, forming a slab inside the box. The system 
is then allowed to evolve for several more orbits, and the 
rate of radial spreading of the slab is used to determine the 
turbulent diffusion coefficient. 

We expect statistical fiuctuations in the evolution of / 
to be large given the chaotic nature of turbulence driven 
by the MRI (Winters, Balbus, & Hawley 2003). Thus, to 
improve our measurement we have performed 22 runs be- 
ginning at different times, reintroducing the contaminant in 
its original distribution each time a run is restarted, and 
evolving the fiow for different periods. We shall group these 
simulations according to the length they were run, label- 
ing each group with different letters, and each run within a 
group by different numbers. Run Al was evolved for one or- 
bit beginning at orbit 16. Runs Bl through B8 were evolved 
sequentially for 0.15 orbits each, starting at orbit 17 (so 
that run Bl provides a measurement of how the contami- 
nant mixes in the fiow during the evolution from 17.00 to 
17.15 orbits; run B2 from 17.15 to 17.30 orbits, etc.). Runs 
CI through C13 were evolved sequentially for 3 orbits each, 
starting at 18.2 orbits (so that run CI is evolved from 18.2 
to 21.2 orbits, etc.). In this way we are able to measure the 
turbulent mixing in the flow over 41.2 orbits of evolution. 

From each run we obtain radial profiles of the az- 
imuthally and vertically averaged contaminant fraction / 
at intervals of 0.05 orbits. For each one of these snapshots, 
wo average the corresponding radial curves over the rmmbcr 
of runs performed, and hence obtain a single curve at each 
time. For example, to obtain the curve that corresponds to 
the radial profile of / at 0.15 orbits, we can average over all 
22 runs (since all were run for at least that long). However, 
for the curve that corresponds to the radial profile of / at 
1.5 orbits, we average over the 13 runs in the group labelled 
C. 



2.3 Calculation of diffusion coefficient 

In order to calculate a value of the diffusion coefficient as- 
sociated with the turbulent flow, we have chosen to solve a 

diffusion equation for the contaminant fraction in the radial 
direction with the same initial condition as that introduced 
in the code: 



dt dx"^ ' 



(7) 



1 \i \x\ < 0.15 
if a; > 0.15 



The solution of this initial value problem is 



f{x,t) = 



1 



'Jx + 0A5^_^Jx-0A5^ 
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(8) 



We then fit this solution to the curves obtained from 
the simulations, treating the diffusion coefficient D as the 
fitting parameter. To test whether D is constant in time. 



the fitting is performed independently for each time in the 
evolution of the contaminant. 

It is worth pointing out that we are using equation (7) as 
a first approximation, and it may not describe the diffusion 
process of the contaminant completely. 



3 RESULTS 

Figure 1 shows the evolution of the / = 0.5 isosurface of 
the contaminant in Run Al at t = 16, 16.05 and 16.1 orbits 
(that is at the start of the run, and 0.05 and 0.1 orbits later). 
At early times, wrinkling of the surfaces on small scales is 
evident, as radial advective transport of the contaminant by 
turbulent eddies occurs. At later times, larger scale distor- 
tions appear, as the contaminant is advected on larger scales 
by larger eddies. 

Figure 2 plots the aaimuthally and vertically averaged 
radial profile of the contaminant from t = 0.05 orbits to 
t = 2.75 orbits after the contaminant is introduced. The av- 
erages are taken from all the data available at each time. 
Over this timescale, steady diffusion of the contaminant to- 
wards a uniform profile (consistent with the radially periodic 
boundary conditions) is evident. 

In Figure 3 we plot fits of the model curve, equation 
(8), to the average radial profiles of contaminant fraction ob- 
tained from the simulations, at the same six times as shown 
in Figure 2. The value of the diffusion coefficient D that pro- 
vides the best fit in each case is shown at the bottom of each 
panel. A total of 57 fits were performed, corresponding to 57 
different instants in the evolution of the contaminant; Fig- 
ure 3 plots six that are representative of the fitting process. 
The error bars shown correspond to the standard deviation 
of the simulation data. 

The values of D obtained by fitting equation (8) ev- 
ery 0.05 orbits of evolution are plotted as a function of the 
time after the contaminant is introduced in the turbulent 
fiow in Figure 4. During the first ~ 0.75 orbits of the evolu- 
tion, D grows from approximately 1.2 x 10~'^i?^/orbit to 
about 5 X 10"^i/^/orbit. After about 0.75 orbits, D be- 
comes more nearly constant; its value does not fall be- 
low « 4 X 10~^ JJ'^/orbit, and does not exceed « 5 x 
lO'^H^/orUt, with an average value between 1 and 3 or- 
bits of (4.7 ± 1.0) X 10"^ii"^ /orbit. As discussed in the next 
section, the change in behavior at 0.75 orbits appears to be 
related to the spatial extent of the contaminant, and the size 
of the largest turbulent eddies in the box. The diffusion time 
tdiff = H^/D, using the saturated value of D at late times, 
is 21 orbits. 

It is useful to compare the magnitude of D measured 
from the simulations with the effective viscosity associated 
with shear stress, using the Prandtl number, PefF = v/D. 
The value of v can be calculated through the relation 



v = acsH = 



« Wx 



(9) 

where cs is the sound speed, << Wxy >> is the time- and 
volume-averaged value of the r-(p component of the stress 
tensor (magnetic plus kinetic), H is the disc scale height, 
and O is the angular frequency. The second equality follows 
from the fact that 

^^ «w » ^^^^ 



Figure 1. Evolution of tiie / = 0.5 isosurface of contaminant in Run A, starting at t = 16 orbits (left panel) when the contaminant is 
introduced into the saturated turbulence, t = 16.05 orbits (center panel), and t = 16.1 orbits (right panel). 



and 

cs = HQ, 



(11) 



Using a value for << Wxy » which is measured di- 
rectly from the simulations (including both the Reynolds 
and Maxwell stress), we obtain v m 0.52i/^/orbit. The 
Prandtl number is then 



Peff = 



v 
D 

0.52 
0.047 
11 



(12) 

(13) 
(14) 



4 DISCUSSION 

Using numerical MHD simulations of turbulence driven by 
the MRI, we have measured the diffusion coefficient for a 
passive contaminant mixed by the flow. When compared 
to the effective viscosity associated with angular momen- 
tum transport, the diffusion coefficient is small, that is the 
Prandtl number Peft = v / D > 1. 

Of course, numerical errors are an important contri- 
bution to the diffusion rate on scales comparable to the 
grid spacing. Ifowcver, on larger (resolved) scales up to the 
box size, there is strong evidence that numerical effects do 
not determine the diffusion rate. The only numerical errors 
in solving equation (6) arise in the advection step of the 



ZEUS code. The magnitude and convergence rate of these 
errors have been studied in Stone & Norman (1992a). For 

the second-order van Leer scheme used in these calculations, 
discontinuities are diffused until they span 10-15 zones in a 
few dynamical times, thereafter very little numerical diffu- 
sion occurs. Figure 2 shows the profiles reach nearly uniform 
profiles in less than 3 orbits, indicative of a much larger diffu- 
sion rate than that associated with numerical errors revealed 
in the advection tests in Stone & Norman (1992a). Thus, 
while numerical diffusion may be important on small scales 
and at early times, it is unlikely to determine the evolu- 
tion observed in Figure 2. Normally the effects of numerical 
diffusion would be assesed by repeating the simulations at 
different resolutions and looking for convergence. However, 
increasing the numerical resolution can also change the prop- 
erties of the turbulence driven by the MRI on small scales; 
such resolution studies are less instructive in this case. 

Figure 4 shows that the best fit value of the diffusion 
coefficient D found by fitting the analytic solution, equation 
(8), to the numerically computed profiles increases nearly 
linearly in time during the first 0.75 orbits of evolution, 
and thereafter saturates at a nearly constant value approxi- 
mately four times the initial value. Thus, the turbulent dif- 
fusion coefficient seems to acquire a meaningful, approxi- 
mately constant value only after the tracer has been spread 
over large scales, comparable to the box size. The possi- 
bility that the diffusion of contaminants in turbulent flows 
follows a non-Fickian diffusion law has been suggested pre- 
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Figure 2. Evolution of the azimuthally and vertically averaged radial profile of contaminant fraction. The corresponding time after the 
introduction of the contaminant into the turbulent flow (in orbits) is shown above each frame. The values of / are shown at each radial 
grid point in the simulation. 



viously by Brandenburg et al. (2004). In our case, the vari- 
ation of D with time may be related to a variation in the 
turbulent velocity fluctuations with scale. Because turbu- 
lence driven by the MRI follows a nearly Kolmogorov scaling 
(HGB), the largest velocity fluctuations are on the largest 
scales. As the profile of the contaminant spreads, it sam- 
ples larger scales, and larger amplitude eddies. Once the 
profile is spread across the largest scale in the simulation 
(the box size), no larger amplitude eddies are possible, and 
D becomes approximately constant. At this stage we have 
not performed simulations with different box sizes; the ef- 
fect of changing the box size may alter the value of D, since 
the amplitude of the MRI turbulence is set by the box size. 
However, the stress would also be modified, and the overall 
efi^ect would be best examined via their ratio. Peg. 

The fact that D varies in time for t < 0.75 orbits calls 
into question the applicability of using equation (7) to de- 
scribe the evolution of the contaminant at early times. A 
more accurate procedure for measuring D would be to use 
the spatial distribution of / at t = 0.75 orbits as the ini- 
tial condition for equation (7), and solve for the resulting 
evolution, so that a solution is found only for times when 
D is known to be constant. Although this may result in a 
slightly modified value for D measured from the simulations. 



our primary result that PeS = v/D > 1 will not change. One 
should also be aware of the possibility that this transport 
process follows anomalous diffusion (Bouchaud & Georges 
1990). 

If the magnetic field acted as a passive scalar in the 
fiow, and one could associate the diffusion coefficient D with 
the coefficient of resistivity r], the magnetic Prandtl number 
would be Pm = v/'q ~ 10, where is the time- and volume- 
averaged sheax stress (Reynolds plus Mcixwell) in the simu- 
lation. In studies of driven turbulence, Yousef et al. (2003) 
found that, for weak magnetic fields, Pm is close to unity. 
There are a number of reasons for this discrepancy. Firstly, 
the magnetic field in shearing box simulations of the MRI is 
clearly not a passive scalar, thus D may be a poor proxy for 
the resistivity. Perhaps just as important, the momentum 
transport rate in the shearing box leads to a much larger 
effective shear viscosity than in driven turbulence simula- 
tions, because the contributions from the Maxwell stress arc 
larger than the Reynolds stress (HGB). Through magnetic 
torques, momentum transport can occur indepent of mass 
motion that would lead to diffusion of passive contaminants. 

The fact that Peff >> 1, so that angular momen- 
tum transport is much more efficient than diffusion of con- 
taminants, has important consequences for the local radial 
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Figure 3. Fits of the theoretical curves [equation (8)] to the numerically measured average radial profiles of the contaminant shown in 
Figure 2. The value of the difi^usion coefficient that provides the best fit is shown at the bottom of each panel, in units of _ff^/orbit. 



transport of small particles embedded in a magnetised disc. 

Global simulations will be required to determine the struc- 
ture and amplitude of radial accretion flows driven by the 
MRI, and therefore assess whether small particles can dif- 
fusive radially outwards in the disc. One of the most im- 
portant questions for global simulations of accretion discs 
is to determine whether large scale poloidal magnetic fields 
can be advected inwards by an accretion flow (leading to 
strong vertical fields in the inner regions of the disc), or 
whether they can diffuse outwards more rapidly than they 
are carried in (Clarke & Pringle 1988). If the diffusion coef- 
flcient D measured here were representative of the vertically 
averaged resistive diffusion coefficient for poloidal magnetic 
fields, then Pm >> 1 would indicate fields would be ad- 
vected inward. However, since poloidal magnetic fields drive 
the MRI, and since reconnection (which is an important sat- 
uration mechanism in the nonlinear regime of the MRI) can 
break fiux-freezing, it is unlikely equation (6) can be used to 
describe the radial evolution of poloidal fields. Global simu- 
lations of thin Keplerian discs which can address these issues 
are needed. 

Acknowledgements: AC acknowledges support from 
CONACYT scholarship 167912, thanks the Department of 
Astrophysical Sciences at Princeton University for their hos- 



pitality, and would like to thank Gordon Ogilvie for fruit- 
ful discussions. ,IS is grateful for financial support from the 
Royal Society, the University of Cambridge, and NSF grant 
AST-0413788 for portions of this work. 



REFERENCES 

Balbus S.A., Hawley J. F. 1991, ApJ, 376,214 

Balbus S.A., Terquem C. 2001, ApJ, 552, 235 

Blaes O.M., Balbus S.A. 1994, ApJ, 421, 163 

Boss A. 1998, Annu. Rev. Earth Planet. Sci., 26, 53 

Bouchaud J. P., Georges A. 1990, Phys. Rep., 195, 127 

Brandenburg A., Kapyla P. J., Mohammed A. 2004, Physics of 

Fluids, 16, 1020 
Clarke C.J., Pringle J.E. 1988, MNRAS, 235, 365 
Fleming T., Stone J.M., Hawley J.F. 2000, ApJ, 530, 464 
Fleming T., Stone J.M., 2003, ApJ, 585, 908 
Gammie C.F. 1996, ApJ, 457, 355 
Glassgold A.E., Najita J., Igea J. 1997, ApJ, 480, 344 
Hawley J. F., Gammie C. F., Balbus S. A. 1995, ApJ, 440, 742 

(HGB) 

Hawley J. F., Stone J.M.. 1998, ApJ, 501, 758 
Heyvaerts J., Priest E.R., Bardou A. 1996, ApJ, 473, 403 
Keller Ch., Gail H. P. 2004, A& A, 415, 1177 
MorfiU G. E. 1983, Icarus, 53, 41 



Diffusion coefficient of a passive contaminant in a local MHD model of a turbulent accretion disc 7 



0,06 



0.05 - 



0,04 



o 



0,02 - 



0.01 



0.00 b 



n 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 — I — r 



<> 

4>-- 



<> 



4> 



t4- 



<> 



<> 



<> 



<>. 



<x><> 



<><> 



<>-- 
--<> 



<><> 



.<x><> 



<><> 



.<> 



<> 



<> 



<> 



<><><> 



<> 



<><> 



<><> 



<> 



<> - 



<> 



<> 



<> 



_] I I L. 



_] I I L. 



_] I I L. 



_] I I L. 



_] I I L. 



_] I I L. 



0,0 



0,5 



1,0 1,5 

Time (orbits) 



2,0 



2,5 



3,0 



Figure 4. Turbulent diffusion coefficient D, obtained from the curve-fitting procedure described in the text, as a function of time. 
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